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1 .0 INTRODUCTION 

A Scout rocket was launched in October 1971 from the National Aeronautics 
and Space Administration (NASA) launching facility at Wallops island, Virginia 
releasing a Bariom ion Cloud (B 1C) over Central America at an altitude of approximately 

20.000 miles. The primary purpose of this experiment was to map the shape of both the 
cloud and the individual magnetic field, lines as the Barium ion Cloud expanded in the 
earth's magnetic field. 

One of the primary tracking .systems for the BIC Project, the Wallops Island 
Image Iniensifier System, was designed and engineered by Electro-optical Systems (EOS), 
a subsidary of the Xerox Company of Pasadena, California. The Image Intensifier System, 
an e iectro-optical device capable of photographing low-light objects by e lecfronica I ly 
enhancing the light received from them, consists of f/l objective lens, an image intensifier 
tube, and a relay lens. Two additional systems, an interference filter for the objective and a 
Flight Research Corporation Mode I 370 Multidata Camera were added at Wallops. 

Si was oniicipaied by WuiJopb personnel that two computer programming systems 
initial ly written by DBA Systems, Inc. for use at the Air Chart and Information Center (ACfC) 
in Saint Louis, Missouri, (1) A Definitive Stellar Camera Calibration program and (2) The 
Geodetic Stellar Camera Orientation program, could be used to calibrate the distortions of the 
Image Intensifier systems and for subsequent data reduction (i.e. determination of direction 
from the tracking station to a specific point on the cloud) of the BIC project photographs. 
However, when stellar calibrations were attempted using the symmetric radial and de- 
centering distortion models, an rms error of photo residuals less than 1 00 micrometers could 
not be achieved. This clearly indicated that the image tube produced a significant level 
of unmodelled systematic error. 

In September of 1971, DBA Systems, Inc. (DBA) received a contract from the 
National Aeronautics and Space Administration (NASA)/Wallops Station, to develop a 
mathematical model to define the Image Intensifier distortions and to implement this 
additional model in the programs identified above on the GE 625 computer at Wallops. 
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In this report we shall outline the development of a suitable mathematical 
procedure for determining the Image Intensifier distortions and the implementation of 
this model in the Wallops computer programs. 
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2.0 ANALYTICAL CALIBRATION OF METRIC CAMERAS 

2.1 Introduction 

A brief summary of the stellar method of calibration of symmetric radial and 
decentering lens distortion simultaneously with elements of interior orientation (x p ,y p/ c) 
and the exterior orientation angles (a, to,*)/ ' s given in this section. 

2 .2 Symmetric Radial Distortion 

The stellar method of simultaneous calibration of symmetric radial distortion 
was developed by (Brown, 1956, 1957, 1964). This work exploited the result from optica! 
ray tracing that the radial distortion 5r of a perfectly centered lens can be expressed as 
an odd powered series of the form: 

6r = K x r 3 + K 3 r 5 + K 3 r 7 +... (1) 

in which, 

i 

r = (x 3 + y 2 ) E = radial distance 

x,y = coordinates of image referred to the principal point as origin. 

Since x,y components of distortion can be expressed as: 

6x = — 6r 

( 2 ) 

Sy = ^- 6r 

r 

it was shown that coefficients of distortion (K % , K 2 , K 3 , . . .) could be introduced directly 
into the projective equations and determined simultaneously with the elements of camera 
orientation in a rigorous least squares adjustment. 
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2.3 


Decentering Distortion 


Only recently has it become appreciated that results for virtually all metric 
cameras are compromised to a small but significant extent by decentering distortion 
(Brown, 1964). Decentering distortion is the result of imperfect centering of lens 
elements. In an earlier form of the Analytical Calibration program, the thin prism 
model was used to account for decentering distortion with satisfactory results. 

According to the thin prism model, a decenfered lens is the equivalent of a perfectly 
centered lens in combination with a thin prism of appropriate deviation and orientation. 

Subsequently (Brown, 1965) it was discovered that a mathematically rigorous 
model for decentering distortion had been developed by Conrady (1919) and that this 
model could be shown to be projective ly equivalent to the thin prism model for first 
order effects but not for higher order effects. Accordingly, an extended form of the 
Conrady model has replaced the thin prism model In the Camera Calibration programs. 

In terms of radial and tangential components, Conrady's model assumes the 

form : 


Ar = 3P r sin (<p- <p 0 ) 
At = P r cos (<p- <p G ) 

in which. 


( 3 ) 


P r = J^r 3 + J 3 r 4 + J 3 r 6 +. . . = profile function of tangential distortion 

< p = angle between positive x axis and radius vector to point x,y 

<p 0 = angle between positive x axis and axis of maximum tangential 
distortion . 
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In Brown (1965) it is shown that Conrady's model can be expressed in terms of x and y 
■ components as: 


Ax = [ P a (r s +2x s ) + 2P 2 xy][l +P 3 r 3 + . . . ] 

Ay = [2P T xy + P 2 (r 3 + 2y 3 )] [1 + P 3 r 3 + ...] 

in which the new coefficients P t , P 3 , P 3 , and P 4 are defined by: 

P i = sin <P 0 

P 2 “ COS 

p 3 = J«A 

?4 = JsA 


(4) 


(5) 


* 


This formulation has the advantage of being a linear expression in the coefficients P^Pg 
when the higher order coefficients P 3 , P 4 are zero. 


2.4 Observational Equations 

The projective equations resulting from an undistorted central projection may be 
written as (Brown, 1957): 


AX + B/i + Cv 
X ~ Xp ~ C 0X+ E n + F v 

A'X + B V + CV 
y ~ y » “ c DX + Efi + Fv 

in which. 


(6) 
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Xp ,y p/C = elements of interior orientation 

\,[l r V = X, Y,Z direction cosines of ray joining corresponding image and 
object points 


~A B C 
A' B' C' 
D E F 


orientation matrix, elements of which are functions of three 
independent angles a,a>,x referred to arbitrary X,Y,Z frame in 
object space. 


If we then let x° , y° represent the observed photo coordinates, the left hand sides 
of the projective equations (6) can be replaced by: 


x-x p - x + v x + x (K x r 3 + K 3 r 4 + K 3 r® + . . .) 

+ [P x (r 2 +2x 3 ) + 2P 3 xy ] [1 + P 3 r 2 + ...] 

vv = v + v + y (K. r 2 + r 4 + K„r 6 + . . .) 

y/p ' y ' % X ^ 

+ [2P a xy + P 2 (r 3 + 2 y 2 )] [1 + P 3 r 3 + .. . ] 


( 7 ) 


in which v and v are photo measuring residuals and, 
k y 



= observed photo coordinates referred to principal point 


r = (x 2 + y 2 ) ~ . 
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3.0 


DEVELOPMENT OF THE IMAGE INTENSIFIER DISTORTION MODEL 


3.1 Introduction 

Development of a mathematical model to define the systematic error of a 
system such as the Image Intensifier requires a comprehensive and uniform set of metric 
data to use as a calibration standard. A stellar photo could provide this standard; however, 
it is almost impossible to acquire an Image Intensifier photograph of a stellar field that 
provides a uniformly dense pattern of stellar images which can be measured with equal 
accuracy. 

Therefore, DBA chose to use a set of photographable targets affixed to an 
ultra-flat 24 inch by 24 inch by 4 inch thick granite surface plate. The relative positions 
of these targets were determined with an accuracy better than 0.0004 inch using DBA's 
proprietary Close Range Analytical Calibration system. A 23 by 23 square array of targets 
spaced at approximately one inch intervals was selected to use as a calibration standard. 

3.2 Calibration Data Acquisition and Red uction 

During the week of 18-22 October 1971, DBA and Wallops Station personnel 
obtained both the close range target calibration and the Image Intensifier (system 
1-09/AC No. 1) photographs of the targeted surface plate. Two sets of acceptable 
Image Intensifier photographs were eventually obtained by stopping the objective to 
f/l 1 and the relay lens to f/4. The first set of three photos were taken at approximately 
2:30 P.M. on 22 October with exposure times of 4,5, and 6 seconds, respectively . Then 
at approximately 6:00P.M. on the same date a final set of 11 exposures were taken at 
15 minute intervals using a 5 second exposure time. There were some 385 target images, 
ranging in quality from poor to good, recorded on the Image Intensifier photographs. 
Typically, the images near the outer edge of the circular format were poor, those in the 
center were fair, while some images in the intermediate area were good. 



Subsequently, these photographs were measured on DBA's 1 micron Mann 
comparator. Then a series of reductions were made in order to determine the random 
errors in the data and the systematic errors due to distortions in both the optical and 
electronic systems of the Image Intensifier,, In this system, of course, the distortions 
of the image tube tremendously outweigh the other two listed sources of error. 

The coordinates of all 529 targets on the surface plate were determined to a 
relative accuracy of better than 0.0004 inch using DBA's Close Range Analytical Cali- 
bration programs. Since this level of error propagates into an error of less than 0 o 75 
micrometers on the photograph, which is an order of magnitude smaller than the 
expected measuring accuracy, the computed coordinates of these targets can be 
considered as being perfect control points for developing an Image Intensifier distortion 
model . 


Four Image Intensifier photographs were selected to be measured and reduced 
for use in model development and evaluation. They consisted of the second frame 
(5~seccnd exposure) from the first set, and frames 2, 6, and 10 of the later set of photos, 
Hereafter, these four frames are referred to as frames 1 , 2, 3 and 4 respectively. 
Primarily frame 1 was used to develop the model and frames 2, 3, and 4 were used to 
evaluate model stability over a long (up to 6 hours) period of continuous operating time. 


These four frames were reduced using DBA's General Multi -frame Analytical 
Calibration program, which (in the single frame mode) employes the identical model as 
that given in Equations (6) and (7), except the direction cosines (A,jz, and v ) are 
computed by: 

. _ XjOG 

* R 




Y-Y c 

R 


v 


Z-Zc 

R 
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X ,Y ,Z 


coordinates of the i ih control target 

coordinates of the center of projection of the Image Intensifier 


in which 


. X c , Y c ., Z c = 


and, 

R = [(X-X c ) 2 + (Y-Y c ) 2 + (Z“Z C ) 5 ]^ . 


This reduction resulted in a set of data (residuals - v„ ,v ) which defined the 
remaining error in all the measurements. Frame 1 was first reduced exercising the symmetric 
radial and decentering distortion models, yielding a residual vector rms of 222.4/im. A 
vector plot of these residuals is shown in Figure 1 . When a similar reduction was made 
without exercising the parameters K 1 ,K 3 ,K g , P a and P 3 of the lens distortion model, a 
residual vector (Figure 2) of 264.5jUm rms was obtained. Subsequent reductions of frames 
2, 3, and 4 (Figures 3, 4, and 5) gave almost identical results with a very similar pattern 
of residual vectors, indicating that, at this level, the distortion errors were essentially 
stable over art extended period of time . 

3.3 The Image Intensifier Distortion Model 

Using the two sets of residual vectors obtained from frame 1, we first applied 
a third degree general polynomial function of the form: 

x - (a o +a^x +a s y +a^xy +a 4 x 3 +a s y 3 
+ a g x y +a 7 xy 3 + a g x 3 +a g y 3 ) 

(8) 

y-(b+bx+by+bxy+bx 2 +by 2 
+ b x 3 y + b 7 xy 2 + b X 3 + b y 3 ) 


A x = 


Ay = 
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In which. 


o Q 

1 

Q 

II 

coefficients of x polynomial 

II 

ja 01 

coefficients of y polynomial 

x,y 

observed image coordinates 

x,y 

true image coordinates (x + v^ , y + v ) 

A*/ ay = 

random plus unmodeled systematic components of v^ 


This procedure met with reasonable success. The rms of residuals (Ax, Ay) were reduced to 
about 50ptm. However, there were significantly large systematic errors remaining, especially 
toward the edge of the photo format. 


Armed with this encouraging result, the general polynomial model was then ex- 
tended to fourth degree (15 terms), and then to a fifth degree function (21 terms) with 
continued success. Surprisingly, the genera! polynomial function consistently fitted the 
residual vectors for the case in which the symmetric radial and decentering distortion 
models were not exercised. With the model extended to a seventh degree genera! polynomial, 
the residual vectors Ax, Ay were reduced to an essentially random pattern with an rms value 
less than 1 0/im. 


At this point, it is necessary to point out some dangers in applying a function 
of this degree. First, unless a large number of image points are used (300 to 400 for 
the initial calibration), the model can become unstable, especial ly with a few "blunder" 
type errors (such as an error in control point identification). Also, with 36 terms being 
determined, the computational effort becomes excessive. . 


Further analysis showed that when the set of data was restricted to image points 
within a circle of radius = 17.5 millimeters from photo center, an extension of Equations (8) 
to fifth degree (using as few as 200 control targets) does an excellent job of removing all 
known systematic errors from the data. Typically, the rms of Ax, Ay residual vectors was 
on the order of 6 to 8 micrometers. 
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FIGURE 1. Plot of residual vectors from frame 1 of targeted surface plate. Symmetric 

radial (K^K^K./ and decentering (P ,P g ) distortion coefficients are applied, 
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FIGURE 3. Plot of residual vectors from frame 2. 
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FIGURE 4. Plot of residual vectors from frame 3. 
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Numerous other mathematical models were investigated but none of them 
approached the effectiveness of the general polynomial function discussed above. The 
model suggested by Wong (1969), for use with television systems, in which a pure 
tangential function was added to the standard optical symmetric radial and decentering 
functions, showed only slight improvement. Interestingly, Wong, Gamble, and Riggins 
(1971) report satisfactory use of a fifth degree polynomial with certain insignificant terms 
eliminated. 

3.4 Evaluation of the Genera! Polynomial Model 

Data from two stellar photographs were used to evaluate effectivensss of the 
general polynomial model. 

First, a stellar photograph from the same Image Intensifier System (1-09/AC No. 1) 
with an exposure time of 0.7 second was selected from a series of stellar photos exposed 
on 18 November 1971 . Approximately 400 stellar images were selected in a nearly 
uniform pattern throughout the photo format. These images were measured and matched 
with SAG stellar catalog coordinates, in preparation for a reduction similar to fhar made 
on the photos from the targeted granite surface plate. The results from this test were 
very similar to those discussed previously. The rms of residual vectors (Figure 6) was 
199.1 p,m. When the seventh degree function was applied the rms (Figure 7) dropped to 
7.4jxm. Application of the fifth degree function to those points within a circle of 17.5 
millimeter radius resulted In an rms (Figure 8) of 7.8/jm. 

Later, when the distortion model was being implemented on the GE-625 computer 
at Wallops Station, a set of actual BIC mission data was obtained to test the program 
modifications. This data was from the Image Intensifier (I-10/H-2) located in Chile 
during the test. The photograph used contained about 350 images, which were mea- 
sured and reduced through stellar identifi cation and updating at Wallops. This reduction 
indicated that system 1-10 has considerably less distortion than system 1-09. The basic 
solution (Figure 9) showed an rms of only 146.7 p, m. Use of the seventh degree function 
dropped the rms (Figure 10) to 6.7 pm, and the fifth degree function limited to points 
within a circle of 17.5 millimeter radius gave an rms (Figure 1 1) of 7.7 micrometers. 
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FIGURE 9. Plot of residual vectors from stellar exposure with Image Intensifier System I- 
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Same as figure 9 after adjustment of points within a circle with 17.5 millimeter 
radius using 5fh„prder general polynomial. 



4.0 


PROGRAM MODIFICATION AND IMPLEMENTATION 


4.1 Introduction 

Two Wallops programs, (1) The Camera Calibration program and (2) The 
Camera Orientation program, required some minor modifications to allow effective 
implementation of the Image Intensifier distortion model. In both programs the changes 
were made such that the program could be used for either a pure optical system such as the 
Wallops Triangulation Camera or the Image Intensifier System. 

Due to basic differences in the logic of the two programs, specifically the 
calibration program uses tape or disk file storage as a data interface between program 
units while the orientation program transfers data through COMMON storage, it was 
necessary to use slightly different versions of the Image Intensifier calibration sub-program. 
Also, the calibration program can use data from an unlimited number of frames to determine 
preliminary distortion coefficients for a single Image Intensifier. 

A description of new control and data parameters, modification to existing 
program units, and implementation of the Image Intensifier distortion program is given in 
the following sections. 

4.2 Calibration Program 

4.2.1 Control Program (MAIN) 

4. 2. 1.1 Program Description 

Program unit MAIN provides sequencing control for the functional units of the 
calibration program. Since significant changes were made in MAIN to allow calibration 
of both optical and Image Intensifier systems, complete listings and flow charts of this 
program unit are given here. 
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Up to 20 systems (this was the existing capability) can be calibrated during 
one computer run. Also, if the Image Intensifler system system is being calibrated, 
data (up to 500 image points) from an unlimited number of photographs (see use of 
control parameters NTYPE and NFRM below) can be used. 

4.2.1 .2 Data 

Card 1 - FORMAT (215) 

NFRM - Number of frames used for current system. NFRM = 1 for 
optical systems 

NTYPE - Type 

0 for optical system 

1 for Image Intensifier 
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4.2.1 .4 Listing 


C CONTROL PROGRAM (MAIN) FOR CALIBRATION 

101 FORMAT (12H1 CALL UPDATE) 

102 FORMAT (12H1CALL AVDIFF) 

103 FORMAT (12H1 CALL PRDISC) 

104 FORMAT (12H1 CALL EMBET ) 

105 FORMAT (12H1 CALL CRVBAL) 

106 FORMAT (12H1 CALL RETURN) 

107 FORMAT (IX 40H 

110 FORMAT (215) 

111 FORMAT (12H1 CALL NCAL ) 

DO 108 1=1,20 
REWIND 7 

READ(5, 1 1 0) NFRM, NTYPE 
DO 115 N=1 , NFRM 
CALL UNIT1 
WRITE(6, 1 01) 

READ (5, 1 07) 

WRITE(6, 1 07) 

CALL UPDATE 
WRITE (6, 1 02) 

CALL AVDIFF 
WR1TE(6, 1 03) 

CALL PRDISC(AAXX) 

115 CONTINUE 

IF(NTYPE .EQ. 0) GO TO 117 
WRITE(6, 111) 

K=-9 
XX = 0. 

WRITE (7) K, XX, XX, XX, XX 
REWIND 7 
CALL NCAL 
GO TO 108 
117 CONTINUE 

IF(AAXX.EQ. 1 .) GOTO 108 
WRITE (6, 1 04) 

CALL EMBET 
WRITE(6, 1 05) 

WRITE (6, 107) 

CALL CRVBAL 
WRITE(6,1 06) 

108 CONTINUE 
STOP 
END 
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4.2.2 


Distortion Calibration Program (PRDiSC) 


4.2.2. 1 Description of Program Changes 

Program unit PRDiSC has been modified to create a data file for input to 
the Image Intensifier calibration program NCAL. Also, the program dimensions and 
loop control parameters have been increased to allow use of as many as 500 control 
points for each frame. Details of these changes are given in the following sections. 

4. 2. 2. 2 Data 

Certain new internal parameters have been added to provide an interface with 
NCAL. These data are described below: 


Name 

Dimension 

Description 

IDXY 

501 

Point identifi cation 

XSV 

YSV 

501 

501 

|$ave X and Y coordinates 

VVX 

VVY 

501 

501 

|$ave X and Y residuals 

NVXY 

NTSV 



Number of points saved 

Number of points saved (temporary) 
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4. 2. 2. 3 List of Program Changes 


Line 

7 

10 

11 

122 

172 

523 

524 

525 

526 

527 

528 

529 

530 

531 

532 

533 

534 

535 

536 

537 

653 

654 

655 

656 

657 

658 

659 


Program Statement 


1 KREJX(501),KREJY(501) 

DIMENSION IDXY(501),XSV(501),VVX(501),VVY(501) 
COMMON/BLKS/XSV, YSV, WX, VVY, IDXY 

DO 44 1=501 

NVXY=0 

NVXY=NVXY+1 
NTSV=NVXY 
IF(SVX*SVY)193,191,193 
191 VVX(NVXY)=0. 

VVY(NVXY)=0. 

XSV (NVXY)=0. 

YSV (NVXY)=0. 

I DXY (NVX Y)=-99 
GO TO 190 
193 CONTINUE 

IDXY(NVXY)=NSTAR 
XSV (NVXY)=TA*1 000. 

YSV (NVXY)=TB*1 000. 

VVX (NVX Y)=VX 
' VVY (NVX Y)=V Y 

DO 705 1=1 , NTSV 
IF(IDXY(I))705,705, 901 

701 WRITE(7)IDXY(I),XSV(I), YSV(I), VVS(I), VVY(I) 

705 CONTINUE 
1=0 

XXX=0. 

WRITE (7)1, XXX, XXX, XXX, XXX 
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4.2.3 


Image Intensifier Distortion Calibration (NCAL) 


4.2.3. 1 Program Description 

Program unit (NCAL) computes the coefficients of an N degree (N<7) general 
polynomial in X and Y. The general form of the equations used; 

X (corrected) = a 0 + a x x + a 3 y + ♦ • .+ a 35 y 7 

Y(correcfed) = b Q + b T x + b 2 y + . . .+ b 35 y 7 

are described in further detail in Section 3.0. 

The program can use data from an unlimited number of frames, however, no 
more than 500 points per frame can be used. Data from each frame is used separately to 
determine a preliminary set of distortion coefficients. Residuals are then computed for 
all points on this frame and compared with a rejection criterion. If a point is either 
rejected or restored the solution is recompute with the adjusted data. When all rejections 
are made for a frame, the normal equation coefficients are accumulated to be used in the 
simultaneous solution using data from all frames. 


4. 2. 3. 2 Data 

Program NCAL control data from three cards and the disk file (7) that was 
formed in program PRD1SC* 

Card 1 - FORMAT (20A4) 

KHEDR - Image Intensifier system identification 


Card 2 - FORMAT 
NPNCH -1° 


(15) 

coefficients are not punched 
coefficients are punched 


Card 3 - FORMAT (15, FI 0.4) 

NX - Number of general polynomial terms to be calibrated (NX- 36) 

i 

RLMT - Points at a distance (X 2 +Y 2 ) 2 > RLMT will not be used in the 
-calibration. 
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Each record on disk file (7) contains the following data: 

a) KPT - point identification; 

b) X - x photo measurement; 

c) Y - y photo measurement; 

d) VX - error in x measurement; 

e) VY ~ error in y measurement. 

A dummy record with KPT=0 terminates the data for each frame and an 
additional record with KPT=-9 terminates the entire set for the current calibration. 
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4. 2. 3. 3 Flow Chari: 


Subroutine NCAL 


START FOR 
DATA FROM 
NEXT PHOTO. 


42-45 


READ DATA 
FROM FILE 7 


KPT(L) -k G, OR 


INITIALIZE 
L = NO. POINTS 
NREJS - 0 
KITER = 0 
WME = 0 



KPT (L) +,0, OR 


COMPUTE 
ERROR VECTOR 
XD AND YD 



m 

























w 

K5 



65 




































4. 2. 3. 4 Listing 


NCAL 

SUBROUTINE NCAL 

COMPUTE IMAGE INTENSIF I ER CALIBRATION COEFFICIENTS 

bth degree general polynomial 

USING MATCHING SETS OF CALIBRATION AND MEAS* COORDINATES 
CALIBRATED COORDINATES in IXCiYC) 

MEASURED COORDINATES IN !X,y) 

INTERNAL COMPUTATION IN M* M* 

DIMENSION KHEDR(aO) iWN( 1E96)/KRJ(B01) 

DIMENSION XC(50l HYC ! 501 H X(5Cl ) / Y< 501 HKpT!501 ) 

DIMENSION SN( 1296),DX(36)/DY(36)^CX(36)/CY{36)/B(36)/WX(36)/WY(36) 

COMMON SN>DX/pY>cXi CY/ WX/WY#WN 
C8MM0N/8LKS/ XC/yC/X/Y/KPTjKRJ 

format statements 

1 FORMAT (1H1) 

2 F 0 RMAT ( 2E 1 6 e 8 ) 

3 FORMAT (515) ^ v ^ 

4 FORMAT (IX, 15/ 2X,E9*4nx>F9*4i2X#F9*4*lX,F9*4* 3X/F9*4i 1X/F9*4# 2X#F6 
#*1/Al/F6*l/Al/3X* F8*l/F£«-1/AHF8*1/Al) 

b FORMAT! Ib/2F 10*4/ 2F 8*1 ) 

6 FORMAT! 1 5 / 2F 1 Q ♦ 4 j 2F 8 * 1 } 

7 FORMAT!/) 

8 FORMAT! 1H1# IXj 5HP0lNT/4X5HX-CALi5X5HY-CALj 7X5HX<-0BS# 5X5HY«*0BS# 

1 6X6hX«C 8MP^ 4X6nY-C0MP>6X2HVX,5X2hVYi 7X5hRD I ST/ 2X6HRAD I AL / 

2 3X5HTANG * / ) 

9 FORMAT ( / 1 8H NO* TERMS USED = 13/ 2QH NO* POINTS USED * 14//) 

10 FORMAT ( 2Ib/5El6«8>2F12»3) 

11 FORMAT ( /22H COMPUTED COEFFICIENTS /) 

12 FORMAT ! 20A4 ) 

14 FORMAT ( / 1 4H MEAN ERROR * F8*i/) 

15 FORM AT ! 2 4H X-DEGREES OF FREEDOM = F6*0/ 24H Y~DEGREES OF FREEDOM 

1= F 6 * C / 24H -X- MEAN ERROR = F8*2/ 24H "Y* MEAN ERR 

2 OR = F6t2//> 

16 FORMAT </22H DEGREES OF FREEDOM » F6,0 /) 

START 

DATA KASTR/KBLNK/1H*/1H / 

LCR = 5 
LPR = 0 
X T 7 = 7 

READ (LCR/12) KHE.DR 
READ ( LCR/ 3 > NPNCH 
READ ( LCR/ 5 ) NX,RLMT 
NXX=NX*NX 
RLMT2 » RLMT*RLMT 
XX a NX 


0 

02 
10 
20 
30 
40 
50 
55 
60 
70 
80 
90 
100 
105 
HO 
120 
130 
140 
150 
160 
i /Q 
180 
190 
200 
210 
220 
230 
240 
250 
260 
270 
280 
290 
30C 
3 1 C 
32C 
33C 
340 
35C 
355 
36C 
37C 
38C 
39C 
40C 
4 1C 
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D0F 3 *(XX*XX) ' 420 

SWME = 0* 430 

SF3 a lOOC* 440 

C SET REF « MATRICES 78 ZERO 45 0 

CALL CLEAR <WN*NXX> 46 ^ 

CALL CLEAR (WX/NX) 47 0 

CALL CLEAR (WY/NX) 48 ^ 

40 CONTINUE 49 ^ 

C READ MEASUREMENT DATA • 500 

06 45 L*l/501 510 

42 READ ( I T 7 ) KPt(L)/ X(L)> Y(L)/ VX/VY 520 

IF l KPT C L ) ) 1 00 j 46/43 530 

43 CONTINUE 540 

XC(L> B X ( L ) - Vx/ SF3 550 

YC(L) = Y(L) • VY/ SF3 560 

KRJ(L) * 1 570 

I F ( X ( L ) * * 2 + Y ( L } **2 ■* RLMT2) 45#45>42 580 

45 CONTINUE 530 

46 CONTINUE 600 

L = L - 1 610 

NREJS =0 620 

KITE* = 0 630 

NME 5 50* 640 

49 CBM I NUE 650 

CALL CLEAR (CX/NX) 660 

CALL CLEAR (CY/NX) 670 

CALL CLEAR. (SN/NXX) 680 

KTEST * 0 690 

KF'LAG * 0 700 

WSQR « (3#0*WME}**2 710 

C COMPUTE NORMAL EQUATIONS AND DEGREES OF FREEDOM 720 

DO 75 K * 1 j L 730 

IF ( KR J ( K ) ) 75/ 75/ 7 1 740 

71 CONTINUE 750 

C COMPUTE B ( I ) MATRIX 760 

B(i) s 1* 770 

B < 2 ) = X ( K ) 780 

B ( 3 ) * Y ( K ) 790 

B { 4 ) a X(K)*Y(K) 800 

B ( 5 ) « X ( K 5 *X ( K ) ‘ 810 

B { 6 ) = Y(K)*Y(K) 820 

B ( 7 ) * B(5)*Y(K) 830 

B(8> = B(6)*X(K) 840 

B { 9 ) a B ( 5 ) *X ( K ) 850 

B(10)= 8 ( 6 ) * Y ( K ) 860 

IF (NX *LE • 10) GO TO 79 870 

B ( 1 1 ) s B ( 9 ) » Y ( K ) 880 

B ( 1 2 ) a B ( 10 ) *X (K ) 890 

8(13)* B(7)*Y(K) 900 
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r 

u 


c 

c 


B ( 14 ) * 

B(9)*X(K) 


B ( 1 5 > = 

B ( 10)*Y(K) 


IF (NX 

*LE • 15) GB 

TO 

B( 16) = 

B( 14 ) *Y(K) 


B( 17)= 

B ( 15 ) *X ( K ) 


B( IS) = 

B(10)*X(K)#X(K) 

B < 19) = 

8 ( 9 ) * Y ( K ) *Y ( K ) 

B ( 20 ) = 

B ( 1 4 ) *X ( K ) 


B ( 21 ) = 

B< 15)*Y(K) 


IF (NX 

«LE< 21) G0 

T0 

B ( 22 ) = 

B ( 20 ) * Y ( K } 


B ( 23 ) ■ 

8< 21 >*X(K) 


B ( 24 ) = 

B(14)*B(6> 


B(25) » 

B ( 5 ) *B ( 15 ) 


B ( 26 ) * 

B ( 9 ) *8 ( 1 0 ) 


B ( 2 7 ) = 

B ( 20 ) * X ( K } 


B ( 2 8 ) = 

B ( 2 1 ) * Y ( K ) 


IF (NX 

• LE * 28) GO 

T d 

B ( 29 ) = 

B ( 27 } *Y ( K ) 


B ( 30 ) = 

B ( 28 ) *X ( K ) 


B ( 3 1 ) » 

B ( 20 ) * B ( 6 ) 


B ( 32 ) = 

8 ( 2 1 ) *B ( 5 } 


B ( 33 ) = 

B ( 1 4 ) * B ( 1 0 ) 


B ( 3 4 ) = 

B ( 9 > *6 ( 15 ) 


8(35) = 

B ( 27 ) * X ( K ) 


B ( 36 ) = 

B ( 28 ) * Y ( K ) 


79 CONTINUE 


XD 

= XC(K) 


YD 

= YC(K) 



call matmfy<b> i>nx, s> i,nx*sn/ 1/ 1 ) 

CALL MATMPYtB* 1/NX/XD# 1/ l* CX* 1*1) 

CALL MATMPY ( B/ 1 * NX/ YD/ 1 1 It CY* 1> 1 ) 

75 CONTINUE 

SET UP XN AND YN AND INVERT 
CALL MATINV(SN/NX>NX/SN) 

CALL MATMPY(SN/NX/NX/CX#NXi \* DX* 0# C) 

CALL MATMPY(SN'NX/NX/CY/NX# 1/DY/C/C) 
CBKPUTE AND PRINT RESIDUALS 

RETURN FQR FINAL RESIDUAL CeMPUTATlQN 
90 CONTINUE 
WVX = 0* 

WVY - 0* 
xx*nx 

DBFX * *XX 
D6FY = «XX 
06 95 K* 1/ L 
B ( 1 ) =1. 

B (2) «? X < K ) 

B ( 3 ) = Y(K) 


910 
920 
930 
940 
950 
960 
970 
980 
990 
1000 
1010 
1020 
1030 
1040 
1050 
1060 
1070 
1080 
1090 
1100 
1110 
1120 
1130 
1140 
1150 
1160 
1170 
1180 
1190 
1200 
1210 
1220 
1230 
1240 
1250 
1260 
1270 
1280 
1290 
1300 
1310 
• 1320 
1330 
1340 
1350 
1360 
1370 
1380 
1390 
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B ( 4 ) =X(K)*Y(K> l^ 00 

B ( 5 ) = X { K ) *X { K ) 14i0 

B ( 6 ) * Y { K ) * Y ( K ) H20 

B { 7 ) = B(5>*Y{K) H30 

B(8) = B ( 6 ) * X { K ) 1^ 4 ° 

B ( 9 ) = B ( 5 ) * X ( K ) 1^50 

B ( 10 ) a B(6)*Y(K) H60 

IF (NX *LE • 10) G0 T0 99 l 470 

B < 1 1 ) = B(9)*Y(K) • .1480 

B(12)= B(10)*X(K) H 90 

B ( 13 ) = B ( 7 ) *■ Y ( K ) 1500 

B( 14)= B(9)*X(K) 1510 

B ( 15 ) s 8(10>*Y(K> 1520 

IF (NX *L E« 15) GO TO 99 i53 ° 

B ( 16 > = B ( 1 4 ) * Y ( K ) 1540 

B ( 1 7 ) « B ( 15 ) *X ( K > I 550 

B(18>= B( 10)*X(K)*X(K) 1560 

B ( 19 ) = B ( 9 ) * Y ( K ) *Y ( K ) 1570 

B ( 20 ) = B(14)*X(K) 1580 

B < 2 1 ) = B ( 1 5 > * Y ( K ) 1590 

IF (NX *LE ♦ 2D G6 TO 99 1500 

B ( 22 ) 3 8 ( 20 ) * Y ( K ) 1610 

B ( 23 ) * B(21>*X(K) 1520 

B ( 24 ) = B<14)*B(6> 1530 

B ( 25 ) = B { 5 ) #8 ( 1 5 } 1640 

B ( 26 ) a B ( 9 ) *B ( 1 0 ) 1650 

B ( 27 ) = B ( 20 ) * X ( K ) 1660 

B(28>» B ( £ 1 ) * Y I K ) 1670 

IF (NX *LE * 28) GO TO 99 1680 

B ( 29 ) = B ( 2 7 } * Y ( K ) 1690 

B(3C)= 8<28)*X(K} 1700 

B(31>» 5 ( 20 ) *B ( 6 ) 1710 

B ( 32 ) = B ( 2 1 ) *B ( 5 ) 1720 

B ( 33 ) = B ( 1 4 ) * B ( 10) 1730 

B ( 34 > = B ( 9 ) #8 ( 15 ) 1740 

B ( 35 ) s B ( 27 ) *X ( K ) 1750 

B ( 36 ) = B ( 2 S ) * Y ( K ) 1760 

99 CONTINUE 1770 

CALL MATMPY(B/1/NX # DXjNX/1#XPI,0/0) 1780 

CALL MATMPY(B/ 1/NXjDY/NX/ l/YPI/C/O) 1790 

VX b(xC(K)-XPI )*SF3 1800 

VY = ( YC ( K ) * YP I )*SF3 1810 

IF(KTEST) 50/50/56 1820 

50 C8NTINUE 1830 

VSGR=VX**2 4 VY**2 1840 

C MAKE REJECT 1 0N OR RESTORATION I860 

IFCKRJ(K) ) 51/51,53 i860 

51 IF(VSGR-InSQR) 52/52/ 56 1870 

C RESTORE 1880 
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52 KRJ(K)*1 
KREJSsNREJS-l 
KFLAG«1 

Ge T8 56 

53 CONTINUE 

IF(VSQR^WSDR) 56i 56/ 54 
C REJECT 

54 KR J { K J s 0 
NREJS^NREJS+I 
KFLAG=1 

56 CONTINUE 

C ACCUMULATE RMS AND DEG* OF FREEDOM 

I F { KR J ( K } ) 62 j 62* 64 
62 CONTINUE 
KX ® KASTR 

KY = KASTR 

GO T9 65 

64 CONTINUE 

kx=kbln K 
ky=kblnk 

DOFX = DOF X + 1* 

DOF Y * D6FY t 1* 

WVX=WVX+VX*VX 

WVY«WVY+VY*VY 

65 CONTINUE 

C PRINT RESIDUALS IF KTEST * 1 

IF(KTEST) 95/95*91 

91 CONTINUE 

RD1ST = SORT ( XP I **2 + YPI**2> 

VRC 3 ( YP I * VY + XPI*VX)/RDIST 
VTC 3 ( XP I *VY - YP 1 * VX ) /RD I ST 

WRITE (LPR/ 4)KP T <K)*XC<K>*YC<K)*XCK>/YCK)*XPI/YPI*VX*KX/VY/KX* 
1 RDIST/VRC/KX* VTC/KX 
IF ( KRj ( K ) ) 95/95*92 

92 CONTINUE 

SWME ■ SWME + VX**2+VY*+2 
DBF 3 DOF + 2* 

XD = XC(K) 

YD 3 YC { K > 

CALL MATMPYCB/ 1/NX*B/ 1/NX/WN/ 1/1) 

CALL MATMPY ( B * 1/NX* XD/ 1/ 1/ W X / 1/ 1 ) 

CALL MATMPY ( 6/ 1/ NX/ YD/ 1/1/ WY/ 1/ 1) 

95 CONTINUE 

WMEX s SGRT< WVX/D6FX) 

WMEY 3 SGRT< WVY/DBFY) 

WME = SGRT ( ( WVX+^VY ) / ( D0FX+DOFY ) ) 

IF (KTEST ) 89/89/88 
88 CONTINUE 

WRITE ( LPR* 7) 


1 89C 
190C 
19 1C 
J92C 
193C 
1 94C 
195C 
196C 
197C 
198C 
199C 
200C 
20 1C 
202C 
203C 
204C 
205C 
206C 
207C 
208C 
209C 
210C 
21 1C 
212C 
213C 
214C 
215C 
216C 
217C 
2 1 SC 
219C 
220C 
221C 
222C 
223C 
224C 
225C 
226C 
227C 
228C 
229C 
230C 
231C 
232C 
2 33C 
234C 
235C 
236C 
237C 


- 38 - 



WR1TECLPR* IE) KHFDR 2380 

WR I TE ( LPR/ 9 ) NX/{_ 2390 

WRITE ( LPR j 7) 2400 

WRITE (LPR/ -15)DGFX/D0FY/WMEX,WhEY 2410 

WR I TE ( LPR/ 14 ) WME 2420 

WRlT£(LPR/7) 2430 

89 CONTINUE 2440 

K 1 TER=K JTER+1 2450 

GO TQ <93/96/96/97/97)/KlTER - 2460 

93 CONTINUE 2470 

G8 T9 49 2480 

96 CONTINUE 2^ 9 Q 

IF(KFLAG) 97/97^49 2500 

97 CONTINUE 2510 

IF(KTeST) 98/98/40 £520 

98 C6NTINUE 2530 

KTEST = 1 2540 

WRITE (LPR/ 8) 2550 

GO T8 90 2560 

C FINAL PROCESSING FOR ALL FRAMES 2570 

100 CONTINUE 2580 

WR I TE ( LPR/ 1) 2590 

WRITEILPR/ 12) KHEPR 2600 

WME * SORT ( SWME/DBF > 2610 

WR I TE { LPR/ 16 ) DQF 2620 

WR I TE { LPR/ 14 ) WMF 2630 

CALL MATINV (WN/NX/NX/SN) 2640 

CALL MAThPY(SN/Nx/NX/WX/NX/ 1/DX/O/C) 2650 

CALL MATMPY(SN/Nx/NX/WY/NX/ 1/OY/O/O) 2660 

WR I TE ( LPR/ 1 1 ) 2670 

NPX=NX*1 2680 

NsQ 2690 

D6 105 M = 1 / NXX/ NpX 2700 

N?N+1 2710 

TA=SGrT(SN(MJ ) * WME/1000* 2720 

TB=DX(N)/TA 2730 

TC=OY(N)/TA 2740 

WR1TE(LPR/ 10) N/M/DX(N>/DY(N)/CX(N) /CY(N) /TA/TB/TC 2750 

IF(NPNCH .EQt 0) GO TO 105 2760 

PUNCH c/DX(N)/DY(M) 2770 

105 CONTINUE 2780 

RETURN 2790 

END 2800 
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4.3 


Orientation Program 


4.3.1 Control Program (MAIN) 

4. 3. 1.1 Program Description 

Also in this set of programs, program unit MAIN serves a control and 
sequencing function when determining the orientation of either an optical camera or 
an Image Intensifier system. The parameter NOJOB continues to control the number of 
jobs that can be processed with a single computer run. However, additional control and 
data parameters, described below, are necessary to implement and refine the (mage 
Intensifier distortion model. 

4.3.1 .2 Data 

Card 1 - FORMAT (NAMELIST/N1/) 

NOJOB - Number of jobs 

Card 2 - FORMAT (15) 

NCAM - Number of systems for which a different set of distortion 
coefficients will be required 

Card 3 - FORMAT (215) 

NTERM - Number of terms in the preliminary Image Intensifier 
distortion calibration 

NTYPF - 0 = optical system 

1 = Image Intensifier 

Card 4 - (3+NTERM) FORMAT (2E16.8) 

XCOEF(I) / NTERM cards containing pre liminary calibration 
YCOEF(!)j coefficients 
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4.3.1 .4 Listing 

C CONTROL PROGRAM (MAIN) FOR ORIENTATION 

COMMON/&LK1/CAM, PL, INST, AZ(2),l:L(2),ROLL(2),ETAC(4, 1),STA 
COMMON/BLK2/ISAT,RMS, ISTR, OBS(3),REF(4),HO, NT, NST 
COMMON/BLK3/IM1 (240), IM2(240), NTC1 (240), NTC2(240), IXA(240), YA 
1 (240) , YA2 (240) , WY(240) , WX (240) , X A2(240) , S L NT (240) 

COMMON/6 LK4/GPHI, GLAM, GAST, STR, KESTH,ESTH,TR, CO, C1,C2,C3 
COMMON/BLK5/GC1 (200), WRA(200), WDE(200),RA(200), DE(200) 
COMMON/6LK6/XPO, YPO,XC, YC, CCPO,SIG(l 0),TK(240), GAM,UWSS 
COMMON/6LK7/AL(2,240), AM(2,240), AN(2, 240),XF(4), YF(4) 
COMMON/BK1/P,T,MF(200),XXPRA, CY,XXPDE,TOBL, FRAC,BDY(16) 
COMMONARUB/LCT,HC,MISS,RANG(240),SLT 
COMMON/T AP/ITAP, STO, IT Y, EQE , LP, NPT (240) 

COMMON/BLKO/NTERM, NTYPE, DXCOEF(36), DYCOEF(36),XCOEF(36), YCOEF(36) 
NAMELIST/N1/NOJOB 
C USE 13 IF ITAP=0 
REWIND 13 
C USE 12 IF IT AP=1 
REWIND 12 

2 FORMAT (2E16.8) 

3 FORMAT (215) 

READ(5,3) NCAM 
DO 500 NC=1,NCAM 
READ(5,3) NTERM, NTYPE 

IF (NTYPE .EQ. 0) GO TO 102 
DO 101 ITN=1, NTERM 

101 READ(5,2) XCOEF(ITN), YCOEF(ITN) 

102 CONTINUE 
ITAP=0 
READ(5,N1) 

DO 300 JOBCT=l, NOJOB 
1 MISS=0 
CALL UNIT1 
CALL UNIT2 



1F(MISS.NE .0) GO TO 300 
CALL UNIT3 

IF (MISS. NE. O.OR. LCT. GE. 15) GO TO 300 
CALL UNIT4 
CALL UNITS 

IF(MISS. NE . 0) GO TO 300 
CALL UNIT7 

. IF(MISS. NE. 0) GO TO 300 
CALL UN1T9 
300 CONTINUE 
500 CONTINUE 
REWIND 13 
REWIND 12 
STOP 
END 



4.3.2 


Subroutine UNIT2 


4.3.2. 1 Description of Changes 

Subroutine UNIT2, Averaging and Differencing of photo measurements, has 
been modified to apply the Image Intensifier distortion corrections to ail stellar and 
object data measurements. These corrections do not apply to the automatic stellar 
search since an alternate set of measured are used by that unit. 


4. 3. 2. 2 Data 

The following data parameters are available in UNIT2 through 
COMM O N/B LK O/ storage . 


NTERM 

NTYPE | 

DXCOEF(36) 
DYCOEF(36) \ 

XCOEF (36) I 
YCOEF(36) \ 


see 4.3.1 .2 


not used here 


see 4.3.1 


O 

• A 
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4. 3. 2. 3 List of Program Changes 


Line Program Statement 

18 COMMON/BLKO/NTERM, NTYPE, DXCOEF(36), DYCOEF(36), 

XCOEF(36), YCOEF(36) 


140 701 IF(NTYPE.EQ.O) GO TO 704 

141 C APPLY IMAGE INTENSIFIER DIST. COR. 

142 AXA=(AXA-XC) *SS2 

143 AYA=(AYA-YC) *SS2 

144 CALL GPCAL(NTERM,XCOEF, YCOEF, AXA, AYA) 

145 AL(1,II)=AXA*SS1+XC 

145 AM(1 , II)=AYA*SS1+YC 

147 GO TO 706 

148 704 AL(1 , ll)=AXA 

149 AM(1,II)=AYA 

150 7 06 CONTINUE 

156 IF(NTYPE.EQ.O) GO TO 6490 

157 C APPLY IMAGE INTENSIFIER CORRECTION 

8 CALL GPCAL(NTERM,XCOEF, YCOEF, AXA, AYA) 

9 AXA=AXA*SS1 

160 AYA=AYA*SS1 

161 GO TO 650 

162 6490 CONTINUE 
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4.3.3 


Subroutine UNIT5 


4.3.3. 1 Description of Changes 


Subroutine UNITS, Preliminary Orientation calibration, has been modified 
to provide the following functions: 

a) apply alternate set of a priori sigmas for principal point (100/^m) for 
Image Intensifier data; 

b) allow the solution to converge (i.e. iterate as many as 8 times if necessary) 
before automatic rejection of point measurements are made; 

c) save x,y photo measurements (X5V,YSV) and error in the measurements 
(VXSV,VYSV) for use in final calibration; 

d) apply final calibration (DXCOEF, DYCOEF) coefficients to ste liar and 
object data measurements. 


4. 3. 3. 2 Data 

The following data parameters are available in Unit 5 through COMMON/ 
BLKO/ and COMMON/BLKS/ storage. 


/BLKO/ 


NTERM 
NTYPE | 

DXCOEF (36) j 
DYCOEF (36) j 

XCOEF(36) 

YCOEF(36) 


see 4.3.1 .2 

final calibration coefficients (for this frame only) 
see 4.3.1 .2 


/BLKS/ 


IDXY (501) 
X5V (501) 
YSV (501) 
VSXV (501) 
VYSV (501) 
KRJXY (501) 


save point identification 
save x measurement 
save y measurement 
save x residual 
save y residual 
•set-up point rejection flag 
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4. 3. 3. 3 List of Program Changes 


Line Program Statement 


13 


COMMON/BLKO/... 

14 


COMMON/BLKS/... 

15 

i 

1 KRJXY(501) 

30 


KFLAG=0 

31 


SS3=1 000. 

32 


SS4=1 001 

33 


W(4, 1)=1 .0D-06 

34 


W(5, 1)=1 .0D-06 

35 


W(6, 1)=1 .0D-06 

36 


XCON=2.5 

37 


IF(NTYPE .EQ. 0) GO TO 503 

38 


XCON=5. 

39 


W(4,l)=l .0D-02 

40 


W(5,l)=l .0D-02 

41 


W (6, 1 )=1 .0D-08 

42 

503 

CONTINUE 

56 


ITE R=0 

57 « 

r* 

START TO ITERATE SOLUTION... 

58 

505 

!TER=!TER+1 

235 


DO 173 M=l,6 

236 


IF(ABS(DD(M,1))-1.0D-07*SQRT(SUMN(M,M)) 

237 

173 

CONTINUE 

238 


GO TO 178 

239 

176 

CONTINUE 

240 


IF (ITER -8) 505,178,178 


178 

CONTINUE 

262 


IF(ABS(UWME-CUWME)-0.5) 5,5,6 

264 


CK=AB$(UWME*XCON) 

280 


IF(NTYPE .EQ. 0) GO TO 315 

281 


IF(KFLAG .NE. 0) GO TO 315 

282 


XPP=XPO*SS2 

283 


YPP=YPO*SS2 

284 


LL=ISTR 

285 


K=LL+1 

286 


DO 301 1=1, K 

287 


KRJXY(I)=1 

288 


IDXY(l)=0 

289 


XSV(l)=0. 
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Line 


Program Statement 


290 

291 

292 

293 

294 

295 

296 

297 

298 

299 

300 

301 

302 

303 

304 

305 

306 

307 

308 

309 

310 

311 

312 

313 

314 C 

315 

316 

317 

318 

319 

320 

321 

322 

323 

324 

325 

326 

327 

328 C 


301 


302 

303 

304 


305 

306 

308 

307 

309 


310 

311 

312 


329 


YSV(1)=0. 

VXSV(!)=0. 

VYSV(i)=0. 

CONTINUE 
DO 305 1=1, LL 
1F(IM1(I)) 306, 305, 302 , 

I F (WX (I) WY(I)) 304,303,304 

KRJXY(l)=0 

IDXY(I)=IM1 (1) 

XSV(I)-(XA(I)-XPP)*SS4 

YSV(I)=(YA(I)-YPP)*S$4 

VXSV(I)=VVX(I) 

VYSV(i)=VVY(l) 

CONTINUE 

!=LL+1 

CONTINUE 

IDXY(l)=-99 

FORMAT (2X, 2I5,2F12.3,2F9.1 , 2F9.4) 

DO 307 K=1 , 1 

WRITE (6, 308) IDXY(K), KRJXY(K),XSV(K), YSV(K), VXSV(K), VY5V(K) 

WX(K),WY(K) 

DO 309 K=1 , ISAT 

WR!TE(6,308) K, IM2(K),XA2(K), YA2(K) 

CALL NCALO 

APPLY ADJUSTED IMAGE INTENSIFIER CORRECTIONS TO STAR MEASUREMENT 
DO 312 K=1 , LL 

1 F (IM 1 (K) .EQ. 0) GO TO 311 

TA=(XA(K)-XPP)*SS1 

TB=(YA(K)-YPP)*SS1 

CALL GPCAL (NTERM, DXCOEF, DYCOEF,TA,TB) 

XA(K)=TA*SS2 
YA(K)=TB*SS2 
IF(KRJXY(K)) 310,310,311 
WX (K)=0. 

WY(K)=0. 

CONTINUE 

WRITE (6, 308) K, IM1 (K),XA(K), YA(K), WX(K), WY(K) 

CONTINUE 

APPLY ADJUSTED IMAGE INTENSIFIER CORRECTIONS TO TARGET 

MEASUREMENTS 

DO 314 K=1 , ISAT 
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Line 

Program Statemenf 

330 

TA=(XA2(K)-XPP)*SS1 

331 

TB=(YA2(K)-YPP)*SS1 

332 

CALL GPCAL (NTERM, DXCOEF, DYCOEF,TA,TB) 

333 

XA2(K)=TA*SS2 

334 

YA2(K)=TB *SS2 

335 

WRITE(6,308) K, IM2(K),XA2(K), YA2(K) 

336 314 

CONTINUE 

337 

XPO=0 

338 

YPO=0 

339 

KFLAG=1 

340 

GO TO 49 

341 315 

CONTINUE 
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4.3.4 


Subroutine GPCAL 


4. 3.4.1 Description 

Subroutine GPCAL is used by subroutine UNIT2and subroutine U NIT5 to apply 
correct both stellar and object data measurements for error due to Image Intensifier 
distortions . 

Using a set of NITERM coefficients, either (XCOEF, YCOEF) or (DXCOEF, 
DYCOEF), and photo measurements (x,y) a corrected set of photo coordinates are 
computed by: 

X (corrected) = a 0 + a 1 x + a 3 y + a 3 xy +...+ a 35 y 7 
Y(corrected) — b 0 + b x x + b s y + b 3 xy +. . .+ b 35 y 7 

where a 0 - a 3S represent DYCOEF or YCOEF. 


Data 


NT 

- number of coefficients in calibration 

DX (36) 

- up to 36 x coefficients 

DY(36) 

up to 36 y coefficients 

XX 

- x photo coordinates (meters) 

YY 

y photo coordinates (meters) 
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4. 3,4. 4 Listing 


SUBROUTINE GPCAL{NT,DX,DY,XX/YY) 

SUBROUTINE T0 APPLY GENERAL POLY* CORRECTIONS - UP TO 5TH DEGREE 
INPUT AND OUTPUT IN METERS * 

INTERNAL COMPUTATIONS IN M • M ? 

DIMENSION B ( 36 ) / DX ( 1),DY( 1) 

SF3= 1CQO ♦ 

X=XX*SF3 
Y = Y Y*SF 3 
B{ 1) =1, 

B(2) =X 
B ( 3 ) => Y 
B(4) = X * Y 
8(5) = X*X 
B ( 6 ) = Y* Y 
B(7) = B(5)*Y 
B(J5) = B < 6 } *X 

er< 9 > = b< 5 )*x ' 

B ( 10 ) = B ( 6 ) ■» Y 

IF (NT .LE* 10) GO TO 99 

B(ll)= B ( 9 ) * Y 

B ( 1 2 ) s B ( 10 ) *X 

B ( 13 ) ■ 3 ( 7 ) * Y 

B< 14) S B( 9)*X 

B( 15)= B( 10)*Y 

IF (NT *LE « 15) g0 TO 99 

B ( 1 6 ) a B ( 1 4 ) * Y 

B ( 1 7 ) * B ( 1 5 ) * X 

B ( 1 8 ) a B ( 1 0 ) *B ( 5 ) 

B ( 1 9 ) = B ( 9 ) *B ( 6 ) 

B ( 20 ) « B ( 1 4 ) *X 

B( 21 ) = B ( 15 ) * Y 

IF (NT *LE* 21) GO TO 99 

8(22)= B ( 20 ) * Y 

6(23)= B<21)*X 

6(24) = B ( 1 4 ) *B ( 6 ) 

B ( 25 ) » B ( 5 ) *B { 15 ) 

B ( 26 ) = B(9)*B(1Q) 

B ( 27 ) = B ( 20 ) * X 

B ( 28 ) = B ( 2 1 ) * Y 

IF (NT -LE* 28) GO TO 99 

B ( 29) = B ( 27 ) *Y 

B { 30 ) = B ( 28 ) * X 

B ( 31 ) = B ( 2 0 ) * B ( 6 ) 

B ( 32 ) = B ( 2 1 ) *B ( 5 ) 

B ( 33 ) = B ( 1 4 ) *B ( 1 0 ) 


1 

c 


1‘ 


e 


i 

c 

ic 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 
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B(3^)= B ( 9 } *B { 15 ) 

B { 35 ) = 8 ( 27 ) *X 
B { 36 ) = 3 ( 28 ) *Y 
99 CONTINUE 

CALL MATMPY(B/1/M,DX/M/1/ X ,0/0) 
CALL MATMPYC8/1/NT/DY/NT/1/ Y ,0/0) 
XX=X/5F3 
YYs Y/SF 3 

return 

END 



4.3.5 


Subroutine NCALO 


4. 3. 5.1 Description 

Subroutine NCALO is identical to program unit NCAL (see Section 4.2.3) 
except (1) the data required in NCALO is stored in COMMON blocks (BLKO) and 
(BLKS) and (2) NCALO uses data from only one photograph to compute a final calibration 
to be used only with that photo. The coefficients obtained here are used by subroutine 
UNIT5, preliminary orientation calibration, to apply the final corrections to both 
stellar and object data (target) coordinates. 


Data 


/BLKO/ 

NX 

- number of terms in calibration 

NTYPE 

not used here 

DX (36) 

final x coefficients 

DY(36) 

- final y coefficients 

DXO(36) 

- not used here 

DYO(36) 

not used here 

/BLKS/ 

KPT (501) 

- point identification 

X (501 ) 

- x photo measurement 

Y(501 ) 

- y photo 

XC(5 01) 

- x (true) 

YC(501) 

- y (true) 

KRJ (501 ) 

- rejection code 
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4.3. 5. 3 Flow Chart 


Subroutine NCALO 


A 


KASTR = * 


INITIALIZE: 


INITIALIZE: 


SET UP ' 

KBLNK = K 


RLMT2 = (RLMT) 2 


REFERENCE MATRICES 



CALIBRATION 

LPR =6 


DOF = -2 (NX) 

S 

WN/WX/WY 


DATA 



SWME = 0 


TO ZERO 























I 


I 




















-z$' 


65 


■ i 
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4.3. 5.4 Listing 


C NCALe c 

SUBR0UT INE NCAL0 0£ 

NCAL F0R ORIENTATION PROGRAM 05 

COMPUTE IMAGE INTENSIFIED .CALIBRATION COEFF IC I ENTS 1C 

5TH DEGREE general POLYNOMIAL 2C 

USING MATCHING SETS OF CAL IBRAT I8N AND MEAS • COORDINATES 3C 

CALIBRATED C0bRD I NATES IN !XC/YC> 4C 

MEASURED COORDINATES IN (X/Y) 5C 

INTERNAL COMPUTATION IN M't M* 55 

6C 

COMMON /BLKB/ NX , NTYPE / DX ( 36 ) / D Y ( 36 ) / DXO < 36 ) / DYO ( 36 ) 7C 

COMMON /BLKB/ KPT<501)^X(501)>Y(501)/XC(50l>/YC(501)>KRJ(50i} 8C 

DIMENSION WN ( 1 296 ) / SN ( 1 296 ) $ KREDR ( £0 ) > CX ( 36 ) / CY ( 36 ) / WX ( 36 ) / 9C 

1 WYI36)#B(36) 95 

COMMON WN/SN/KHEdR/CX/CYjWX/WY/XPI/ YPI IOC 

C FORMAT STATEMENTS 1 1C 

1 FORMAT ( 1H1 ) 1 2C 

2 F 6 RMAT ( 2E 1 6 • 8 ) 13C 

3 FORM Af ( 5 I 5 ) 14C 

4 FORMAy! IX/ I5/2X/F9t4/ 1X,F9*4/2X/F9.4* 1X/F9.4/3X/F9.4/ 1X/F9.4/ 2X/F6 15C 

*tl/Ai/F6il/Al/3X/ F8«l/F&*l/Al/F8*l/Al) 16C 

5 F0RMATU5/2F1O*4,2F8»1 ) 17C 

6 FORMAT! 15, 2F10* 4/ 2F8 > 1 ) ISC 

7 FORMAT!/) 19C 

8 FORMAT ! iHl/ 1X/5HPBINT/4X5HX-CAL/5X5 HY*CALj 7X5HX-0BS, 5X5HY-0BS, 20C 

1 6X6HX-C0MP/ 4X6HY-C0MP/6X2HVX/5X2HVY/7X5HRDIST/2X6HRADIAL/ 21 C 

2 3X5HTANG./) 22C 

9 FORMAT ! / 18H NO • TERMS USED = 13/ 2QH NO. POINTS USED = 14//) 23C 

10 F0RMAK2I5/5E16.8/2F12.3) 24C 

11 FORMAT ! /22H COMPUTED COEFFICIENTS /) 25C 

12 FORMAT ( 20A4 ) 26C 

14 FORMAT (/ 1 4H MEAN ERROR « F8*l/) 27C 

15 FORMAT ! 24H X-DEGREES OF FREEDOM = F6«0/ P4H Y-DEGREES OF FREEDOM 28C 

1* F 6 • 0 / 24H -X- MEAN ERROR ? F8*2/ 24H MEAN ERR 29C 

20R » F8.2//> 30C 

16 FORMAT ! /22H DEGREES OF FREEDOM = F6.0 /) 31C 

C START ‘ 32C 

DATA KASTR/KBLNK/1H*/1H / 330 

LCR=5 - 340 

LPR«6 35C 

RLMT = 25 * 360 

NXX = NX*NX 390 

RLMT2 = RLMT*RLMT 400 

XX s NX 410 

DOF » -(XX+XX) 

SNME « 0* 430 
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I 


SF3*10CC. 44C 

C SET REF. MATRICES T9 ZERO 45C 

CALL CLEAR (WN/NXX) 46C 

CALL CLEAR (WX/NX) ' 470 

CALL CLEAR (WY/NX) 480 

40 CSNTINLE 490 

C SET-UP INPUT DATA 500 

06 45 L«l/501 510 

I F { KPT ( L ) ) 46/44, 43 530 

43 CONTINUE 540 

XC(L> - X(L) + XC<L)/SF3 550 

YC { L ) = V{L) + YCID/SF3 560 

KRJ { L } « 1 570 

IF(X<L>**2 + Y ( L ) **2 - RLMT2 ) 45/45,44 580 

44 KRJ(L) = 0 585 

45 CONTINUE 590 

46 CONTINUE 600 

IF C L c LE « NX) G6 TO 200 605 

L = L-l 610 

NREJS = 0 620 

KITER = 0 630 

WME - 50* 640 

4? CONTINUE 65Q 

CALL CLEAR (CX/NX) 660 

CALL CLEAR (CY*NX> 670 

CALL CLEAR(5N/NXX) 680 

KTEST = 0 - 690 

KFLA6 = 0 700 

VJSQR = ( 3 « 0*WME ) **2 710 

C C8MPUTE NORMAL EgUaTIQNS AND DEGREES 0F FrEEDBM 720 

D6 75 K= 1/L 730 

IF ( KR J ( K ) ) 75/75/71 740 

71 CONTINUE 750 

C COMPUTE B(I) MATRIX 760 

B { 1 ) « 1. 770 

B { 2 ) a XU) 780 

B C 3 ) = Y ( K ) 790 

B { 4 ) = XU)*Y(K) 800 

B ( 5 ) a X ( K ) *X ( K ) 810 

B ( 6 ) » Y U ) * Y ( K ) 820 

B < 7 ) = B(5)*Y(K) 830 

B { 8 ) a B(6>*X(K> 840 

B(9) a B ( 5 ) * X ( K ) 850 

B ( 10 ) » B ( 6 } * Y { K ) • 860 

IF (NX *LE • 10) GO T9 79 870 

B { 1 1 ) s B ( 9 ) * Y { K ) 880 

B( 12)= B( 10)*XU) 890 

B ( 13 ) = B ( 7 ) * Y ( K ) 900 

B ( 1 4 ) ® B{9)*X(K) 910 
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c 


c 

c 


79 


75 


90 


B { 15 ) = B(10)*Y<K) 

IF (NX *LE * 15) 68 TB 79 
B < 16)= B ( 1 4 ) * Y ( K ) 

B( 17)= B ( 1 5 ) *X ( K ) 

B ( 1 8 ) = B ( 10 ) *X ( K ) *X ( K ) 

B ( 19 ) = B(9)*Y(K)#Y(K) 

B ( 20 ) = B(14)*X(K) 

B ( 2 1 ) = B ( 1 5 ) * Y ( K ) 

IF (NX *LE> 21) G0 TB 79 
B ( 22 ) = B { 20 ) * Y ( K ) 

B ( 23 ) = B ( 21 > *X ( K } 

B ( 24 ) = B ( 14 ) *B ( 6 ) 

B ( 25 ) = 8 ( 5 ) *B ( 15 ) 

B ( 26 ) = B ( 9 ) *B ( 1 0 ) 

B ( 27 ) = 8 ( 20 ) * X ( K ) 

B ( 28 ) = B ( 2 1 ) * Y ( K ) 

IF (NX «LE* 28) 68 T0 79 
B ( 29 ) = B ( 27 ) * Y ( K ) 

B ( 30 ) = B ( 28 ) *X ( K 5 
B(3i)= 8 ( 20 ) *8 ( 6 ) 

B ( 32 ) = B ( 21 > *B < 5 > 

B ( 33 ) = B(14)*0(lO) 

B ( 35 ) » B(27)*X(K) 

B ( 36 ) = B(28)*Y<K) 

CONTINUE 
XD * XC ( K ) 

YD * YC(K) 

CALL MATMPYtBi l/NX,Bi 1 jNX/SN> 1) 

CALL MATMPY(B/1jNX#XD#1j1jCX# 1# 1) 

CALL MATNPY ( B/ X, nX> YD/ 1 / 1 > C Y> 1 - 1 ) 
CONTINUE 

SET Up XN AND YN AND INVERT 
CALL NATINV(SN/NX/NXjSN) 

CALL MATMPY(SMNX/NXiCX/NX/ 1/DX/O/C) 

CALL MATMPY(SN/NXiNX/CYjNX/ 1/DY/C/C) 
compute AND PRINT RESIDUALS 

return for FINAL RESIDUAL COMPUTATION 
CONTINUE 
WVX =« 0* 

WVY =? C« 

XX = NX 

D8FX e *XX 
D8FY = *XX 
D8 95 K= 1/L 
B ( 1 ) = 1 . 

B ( 2 ) = X ( K ) 

B ( 3 ) =Y(K) 

B ( 4 ) =X ( K ) *Y ( K ) 


92C 
93C 
94C 
95C 
96C 
97C 
98C 
99C 
100C 
101C 
102C 
JC3C 
104C 
1 Q5C 
106C 
107C 
10 SC 
1C9C 
HOC 

me 

112 C 

113C 

i \ ii.r 

115C 

116C 

U7C 

1 1 sc 

119C 
120C 
121C 
1 22C 
1 23C 
124 C 
125C 
126C 
127C 
128C 
1H9C 
130C 
1 3 1C 
1 32C 
1 33C 
1 34C 
135C 
136C 
1 3 7 C 
138C 
139C 
140C 
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B ( 5 ) =. X ( K } * X ( K ) 141< 

B ( 6 ) = Y ( K ) * Y < K ) 142f 

B ( 7 ) a B ( 5 ) *Y ( K ) 143C 

B ( 8 ) = B(6)*X(K) 144( 

B ( 9 ) = B(5)*X(K) 145C 

B ( 1 0 ) = B(6)*Y(K) i46f 

IF (NX *LE* 10) G0 TO 99 147C 

B ( 1 1 ) a B(9)#Y(K) 1 48C 

B ( 12 ) a B( 10)*X<K) 1 49C 

8(13)= B ( 7 ) *Y ( K ) 150C 

8(14)= B(9)*X(K) 1510 

B ( l5 > c B(10)*Y<K) 152C 

IF (NX *LE • 15) T9 99 153C 

B ( 16) = B ( 14 ) * Y ( K } 1 54C 

8(17)= 8 ( 15 ) *X ( K ) 155C 

B ( 18 ) = B( 10 ) * X ( K ) *X ( K ) i 560 

B ( 19 ) = B ( 9 ) * Y ( K ) * Y ( K ) 1 57C 

B ( 20 ) = B(14)*X(K) 158C 

B ( 2 1 ) = 8 ( 15 ) *Y ( K ) 159C 

IF (NX *LE * 21) 66 T9 99 160C 

B ( £2 ) s B ( 2 0 ) * Y (' K ) 1 6 1 C 

B ( 23 ) = 8(21 ) * X ( K ) 1 62C 

8(24)= B ( 1 4 ) * B ( 6 ) a 630 

B(25)= 8<5>*5(15) 164C 

8(26)= B ( 9 ) #B ( 1 0 ) 1 65 C 

B ( 27 ) = B(20)*X(K) 166C 

B ( 28 ) = B(21)*Y(K) 167C 

IF (NX 28 ) 68 T8 99 - 168C 

B ( 29 ) = 8 ( 27 ) * Y ( K ) 169C 

6(30)= 8 ( 28 ) *X ( K ) 170C 

B ( 3 1 ) = B ( 20 ) *B ( 6 ) 1 7 1 C 

B ( 32 ) = B < 2 1 ) *B ( 5 ) 172C 

B ( 33 ) = B ( 1 4 } *B ( 10) 1730 

8(34)= 8 { 9 ) *B ( 1 5 ) 174C 

B ( 35 ) = B ( 27 } *X ( K ) 175C 

B ( 36 ) = 9(28)*Y(K) 176C 

99 CONTINUE 177C 

CALL MATMPY(B/l/NX,DX/NX/l/XPI,0/0> 178C 

CALL MATMPY(B/l/NX,DY/NX/l/YPI/0/0) 179C 

VX =(XC(K)-XPI ) *SF3 180C 

VY = ( YC ( K ) * YP I )*SF3 1 S 1C 

IF(KTEST) 50/ 50/56 182C 

50 CONTINUE 183C 

VSQR 5 VX**2 + VY**2 184C 

C MAKE REJECTION 6R RESTORATION 185C 

IF<KRJ(K>) 51/51,53 186C 

51 IF(VSGR-WSCR) 52,52/56 1S7C 

C RESTORE 1S8C 

52 Krj<K)=i 18SC 
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c 


c 


c 


NRE JS=NRE JS-1 
KFLAG=1 
GQ T9 56 

53 CQNTINUE 

IF { VSQFc^wSQR ) 56/56/54 
REJECT 

54 KR J ( K ) = 0 
NREJS=KREJS+1 

kflag= 1 

56 C0NTINUE 

ACCUMULATE RMS AND DEG* QF FREEDOM 
IF(KRJ(K) ) 62/62/64 
62 CONTINUE 
KX * KASTR 
Ky b KASTR 
G0 T9 65 

64 CONTINUE 
KXtKBLNK 

ky=kbl nk 

D9FX = D9FX + 1* 

D8FY = DOF Y + 1# 

wvx»wvx+vx*vx 

WVY«WVY+VY*VY 

65 C0NTINUE 

PRINT RESIDUALS IF KTEST = 1 
IF(KTEST) 95/95/gi 

91 CQNTINUE ’ 

ROIST s SORT ( XP I * * 2 + YPI**2) 

VRC 3 { YP I * VY + XP 1 * VX ) /RD I ST 
VTC 3 ( XP I * VY - YPI*VX)/RDIST 

WRITE (LPR/ 4)KPt(K)/XC(K)/YC(K)/X(K)/Y(K)/XPI/YPI/VX/KX/VY/KX/ 
1 RDIST/VRC/KX/ VTC/KX 
J F ( KR J { K ) J 95/95/92 

92 CONTINUE 

SWME = SWME + VX**2+VY**2 
DQF 3 D9F + 2* 

XD 3 XCfK) 

YD = YC<K) 

CALL MATMPYIB/ l/NX/8/ 1/NX/WN/ 1/ 1 ) 

CALL MATMPY ( B/ 1 / NX/ XD/ 1 / 1 / WX/ 1 / 1 ) 

CALL MATMPYIB / 1/NX/YDi i> 1/WY/ 1/ 1 ) 

95 CQNTINUE 

WMEX 3 SQRT(WVX/06FX) 

WMEY 3 SQRT( WVY/08FY) 

WME 3 SGRT < ( WVX+wVY ) / ( D6FX+D6FY) ) 

IF(KTEST) 89/89/88 
88 CQNTINUE 

WRITE (LPR/ 7) 

WR I TE { LPR/ 12 ) KHEDR 


1900 

1910 

1920 

1930 

1940 

1950 

1960 

1970 

1980 

1990 

8000 

2010 

2020 

2030 

2040 

2050 

2060 

2070 

2080 

2090 

2100 

2110 

2120 

2130 

2140 

2150 

2160 

2170 

2180 

2190 

2200 

2210 

2220 

8230 

2240 

225C 

2260 

2270 

2280 

2290 

2300 

2310 

2320 

2330 

2340 

2350 

2360 

2370 

2380 
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WR 1 T£ { LPR# 9 ) NX/l 
WRITE (LPR# 7) 

WRITE (LPR# 15)D6FX#D0FY#WhEX#wMEY 
WR I TE ( LPR# 14 ) WME 
WR I TE ( LPR/ 7 ) 

89 CONTINUE 

KITERsKITER+1 

Ci© TQ <33#96#96#97#97)#KITER 
93 CONTINUE 
GO T© 49 

96 CONTINUE 
IF(KFlAG) 97#97#/*9 

97 CONTINUE 
IF(KTEST) 98/98>l00 

98 CONTINUE 
KTEST = 1 
WRITE ( LPR# 8) 

GO T0 90 

C FINAL PROCESSING FOR ALL FRAMES 

100 CONTINUE 

WRITEILPR/ 1) ■ 

WR I TE ( LPR# 12 ) KHEDR 
WME = SORT < SWME/DOF ) 

WRITE<LPRil6> 08F 
WRITER LPR# 14) WME 
CALL MATlNV(WN/NXiNX#SN> 

CALL MATMPY(SN/NX#NX#WX#NX# 1/DX#0#C) 

CALL MAThPY(SN/NX#NX#WY#NX# 1#DY#0#C) 

WRITEILPR# 11 ) 

NPX=Nx+l 
N= Q 

DO 105 M=1#NXX#NPX 
N« N+ 1 

TA=SQRT(SN(M) ) * WME/1000. 

T8=DX(N)/TA 

TCsDY{N)/TA 

WRITE (LPR# 10) N,m#DX (N)/DY(N)/CX(N) /CY(N) #TA/TB#TC 
105 CONTINUE 
200 C9NTINUE 
RETURN 
END 


239C 
240C 
24 1C 
242C 
243C 
244C 
245C 
246C 
247C 
248C 
249C 
250C 
25 1 C 
252C 
253C 
254C 
255C 
256C 
257C 
258C 
259*' 
260( 
26 1C 
262 f 
263( 
264C 
265i 
266( 
267c 
268( 
269c 
270c 
27 1< 
272 < 
273c 
274c 
275c 
276c 
277c 
279 [ 
280c 
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4.4 


General Subroutines 


4.4.1 Introduction 

Three general subroutines CLEAR, MATMPY, andMATINV were used in this 
set of programs. The function and listings of these subroutines are included for reference 
only . 

4.4.2 Subroutine CLEAR 

Subroutine CLEAR set (N) elements of a real array (X) to zero. 


C 

SUBROUTINE CLEAR { X* N ) 

c clear storage re zero 

D I KENS I BN X ( 1 > 

DO 1 IrliN 
1 X ( I ) = 0 • 

RETURN 

END 


2810 

2820 

2830 

2840 

2850 

2860 

2870 

2880 
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4o4.3 Subroutine MATMPY 


Subroutine MATMPY multiplies a matrix B with dimensions of NRB rows and 
NCB columns by matrix A with dimensions of NRA rows and NCA columns and stores the 
results in Co Codes for optional multiplications and accumulations are given in the listing. 


C 


SUBR9UT INE MATMR V < A# NR A* NC A; B* NRB* NCB* C* Ml > M2 ) 


MATMPY A FORTRAN MATRIX MULTIPLY (AND ADD) PACKAGE 

will Handle arrays of any dimensions 

NO CHECK IS MADE TO ENSURE A VALID PRODUCT OR SUM* 
ES FOR OPTIONAL MULTIPLIES AND ADD I T 1 &NS 


C 


10 


20 


30 


40 


50 


Ml 


R 

0 


A*B 

1 


AT*B 

2 


A*BT 

3 


AT*BT 

DIMENSION 

A ( ; 

IF (Ml 

-1) 40 

*50*10 

NCC 5 

NRB 


I NR c 

NCB 


MB = 

1 


I NCB 

« NRB 


IF(M1 

-2) 50 

*30*20 

NRC = 

NCA 


INCA 

8 1 


MA = 

NRA 


00 T8 

60 


NRC * 

NRA 


INCA 

8 NRA 


MA ■ 

1 


GO T9 

60 


NRC 5 

NRA 


NCC =? 

NCB 


MA 8 

1 


MB s 

NRB 


INR 3 

NCA 


INCA 

8 NRA 


INC8 

8 1 


GO T0 

60 


NRC 3 

NCA 


NCC 3 

NCB 


INR = 

NRA 


MA 8 

NRA 


MB 3 

NRB 


INCA 

8 1 


INCB 

8 1 



ME 

0 

1 

1 

*2 


C 

R 

C + R 
C " R 


-R 


b ( i ) / cm 


289C 
290C 
291C 
292C 
293C 
294C 
295C 
296C 
297C 
298C 
299C 
300C 
30 1C 
302C 
303C 
304C 
305C 
306C 
307C 
308C 
309C 
3 1 0C 
3 1 1C 
312f 
313C 
314( 
315C 
3 1 6 f 
317c 
318c 
319c 
320 f 
321' 
322' 
323' 
324i 
325' 
326' 
327' 
328' 
329' 
330' 
331' 
332' 
333- 
334' 
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POSITIVE NEGATIVE 

3350 


MR RP=R RP=*R 

3360 


MC C = C + RP C = RP 

3370 


A 

3380 

60 

IF ( M2 ) 70# 80# 90 

3390 

70 

IF ( M2+ 1 ) 110# 100/80 

3400 

80 

MO - 1 

3410 


MR « 1 

3420 


GO T0 120 

3430 

90 

MC 3 1 

3440 


MR 3 1 

3450 


68 TO 120 

3460 

100 

MC 3 1 

3470 


MR a -1 

3480 


68 T8 120 

3490 

110 

MC * ^1 

3500 


MR = -1 

3510 

120 

CONTINUE 

3520 


D6 190 I 3 1 # NRC 

3530 


I J 3 I 

3540 


INTLA * ( I-1)*MA + 1 

3550 


D6 180 J® 1# NCC 

3560 


D A - 

\ \ w " 

3570 


IB = <J^1)*MB + 1 

3580 


IA 3 INTLA 

3590 


D6 130 K= 1# INP 

3600 


R = R + A( I A) «B ( IB ) 

3610 


IA 3 I A + INCA 

3620 


IB * IB + I NCB 

3630 

130 

CONTINUE 

3640 


IF { MR ) 14Q# 140/ 150 

3650 

140 

Rs -R 

3660 

150 

IF(MC) 160# 160# 170 

3670 

160 

C ( I J ) = R 

3680 


G0 T9 180 

3690 

170 

C(IJ) = C(I J) + R 

3700 

180 

IJ 3 IJ + NRC 

3710 

190 

C6NTINUE 

3720 


return 

3730 


END 

3740 
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4.4.4 


Subroutine MATINV 


Subroutine MATINV inverts a matrix with dimensions of NROW rows and 
NCOL columns and stores the inverse matrix in B. A and B may be the same array. 


C 3750 

SUBROUTINE MATINV ( A/NR6k, NC8L* B ) 3760 

DIMENSION a<1)/B<1) 3770 

Nr=NR0W 3780 

NC=NC0L 3790 

NA=NR*NC • 3800 

DG 5 I=1/NA 3810 

5 B ( I ) = A ( I ) 3820 

DG 25 J=l, NR 3830 

I + NR* { I 1 ) 3840 

C*i*C/B(N) 3850 

B(N)*1*C 3860 

DG 10 N=I,NA/NR 3870 

10 B ( N ) s C * B ( N ) 3880 

DG 25 K = 1 i NR 3890 

N^K+NR* ( I i ) , 3900 

IF ( I -K ) 1 5 j 25/ 15 3910 

15 C c B { N ) 3920 

B ( N ) =0 * C 3930 

DG 20 U=i/NC 334 0 

l?NR*<J-l) 3550 

N = K + L 3960 

l »I+U 3970 

20 B { N ) =B ( N ) *C*B ( L ) 3980 

25 CONTINUE 3990 

RETURN 4C00 

END *010 
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